subroutine calculate_characteristics_anomal(points,time,edges,nodesforedges,&
    vertex3D,nVer3D,numbersEl3D,igT,jgT,ggT,nt,nodesCount3D,elemsCount3D,&
    edgesCount3D,edgesTermCount3D,edgesNotTermCount3D,mTerm3d)
    common /characteristics/ nPoints
    real*8 time,points,strength,coef,vertex3D,aX,aY,x,y,z,ggT
    integer timeLevel,edges,elemsCount3D,edgesCount3D,searchElement3D,edgesTermCount3D,&
    edgesNotTermCount3D
    character*4 Ext
    dimension time(nt),points(3,nPoints),strength(3),coef(3),edges(12,elemsCount3D),&
    nodesforedges(3,edgesCount3D),vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),&
    aX(3),aY(3),numbersEl3D(12,nodesCount3D),igT(edgesTermCount3D+1),jgT(mTerm3d),&
    ggT(mTerm3d)

    if(nPoints.ne.0) then
        call makedirqq('.\results\E')
        do j=1,nPoints
            x = points(1,j)
            y = points(2,j)
            z = 0.d0
            nEl = searchElement3D(vertex3D,edges,nodesforedges,x,y,z,nodesCount3D,&
            elemsCount3D,edgesCount3D)
            write(Ext,'(i4.4)') j
            open(500,file='.\results\E\Eax.'//Ext,mode='write')
            write(500,"(/,A)") '      t (mc)             ExA (mB)'
            open(501,file='.\results\E\Eay.'//Ext,mode='write')
            write(501,"(/,A)") '      t (mc)         ExA (mB)'

            do timeLevel=10,nt
                call GetTimeCoef(time,coef,nt,timeLevel)
                call GetSolutinonToNode(vertex3D,nVer3D,edges,nodesforedges,numbersEl3D,&
                Ax,Ay,nEl,nodesCount3D,elemsCount3D,edgesCount3D,timeLevel,igT,jgT,ggT,&
                edgesTermCount3D,edgesNotTermCount3D,mTerm3d,x,y,z)
                call GetStrengthA(aX,aY,strength,coef)
                write(500,*)time(timeLevel)*1000.d0,strength(1)*1000.d0
                write(501,*)time(timeLevel)*1000.d0,strength(2)*1000.d0
            enddo
            close(500)
            close(501)
        enddo
    endif
end
!*************************************************************************
subroutine GetTimeCoef(time,coef,nt,timeLevel)
    dimension time(nt),coef(3)
    real*8 time,coef,t1,t2,t3,dt,dt1,dt0
    integer timeLevel

    t1=time(timeLevel-2)
    t2=time(timeLevel-1)
    t3=time(timeLevel)
    dt=t3-t1
    dt1=t2-t1
    dt0=t3-t2
    coef(1)=dt0/dt/dt1
    coef(2)=-dt/dt1/dt0
    coef(3)=(dt+dt0)/dt/dt0
end
!*************************************************************************
subroutine GetSolutinonToNode(vertex3D,nVer3D,edges,nodesforedges,numbersEl3D,&
    aX,aY,nEl,nodesCount3D,elemsCount3D,edgesCount3D,timeLevel,igT,jgT,ggT,&
    edgesTermCount3D,edgesNotTermCount3D,mTerm3d,x,y,z)
    real*8 vertex3D,a,ggT,getValueOnEdge,aX,aY,temp,x,y,z,hx,hy,hz,x1,x2,&
    y1,y2,z1,z2,newX,newY,newZ
    integer elemsCount3D,edgesCount3D,edges,timeLevel,edgesTermCount3D,&
    edgesNotTermCount3D
    character*4 Ext
    dimension edges(12,elemsCount3D),nodesforedges(3,edgesCount3D),&
    vertex3D(3,nodesCount3D),nVer3D(14,elemsCount3D),aX(3),aY(3),&
    numbersEl3D(8,nodesCount3D),igT(edgesTermCount3D+1),jgT(mTerm3d),&
    ggT(mTerm3d),a(12,3),temp(3)

    nFile = 502
    do i=1,3
        write(Ext,'(i4.4)') timeLevel+i-4
        open(nFile,file='.\calculations\a\A.'//Ext,access='direct',recl=8,mode='read')
        do j=1,12
            a(j,i) = getValueOnEdge(igT,jgT,ggT,edges(j,nEl),edgesTermCount3D,edgesNotTermCount3D,mTerm3d,nFile)
        enddo
        close(nFile)
        aX(i)=0.d0
        aY(i)=0.d0
    enddo

    hx = vertex3D(1,nVer3D(2,nEl))-vertex3D(1,nVer3D(1,nEl))
    hy = vertex3D(2,nVer3D(3,nEl))-vertex3D(2,nVer3D(1,nEl))
    hz = vertex3D(3,nVer3D(5,nEl))-vertex3D(3,nVer3D(1,nEl))
    x1=vertex3D(1,nVer3D(1,nEl))
    x2=vertex3D(1,nVer3D(2,nEl))
    y1=vertex3D(2,nVer3D(1,nEl))
    y2=vertex3D(2,nVer3D(3,nEl))
    z1=vertex3D(3,nVer3D(1,nEl))
    z2=vertex3D(3,nVer3D(5,nEl))
    call transformationCoordinats(x,y,z,hx,hy,hz,x1,x2,y1,y2,z1,z2,newX,newY,newZ)

    do k=1,3
        do i=1,12
            call  basisFunctionOnVectPar(i,newX,newY,newZ,hx,hy,hz,a(i,k),temp)
            aX(k) = aX(k) + temp(1)
            aY(k) = aY(k) + temp(2)
        enddo
    enddo
end
!*************************************************************************
subroutine GetStrengthA(aX,aY,strength,coef)
    real*8 aX,aY,strength,coef
    dimension aX(3),aY(3),coef(3),strength(3)
    
    strength(1)=0.d0
    strength(2)=0.d0
    do i=1,3
        strength(1)=strength(1)-coef(i)*aX(i)
        strength(2)=strength(2)-coef(i)*aY(i)
    enddo
end
!*************************************************************************
real*8 function getValueOnEdge(igT,jgT,ggT,node,edgesTermCount3D,edgesNotTermCount3D,mTerm3d,nFile)
    real*8 ggT,v
    integer edgesTermCount3D,edgesNotTermCount3D
    dimension igT(edgesTermCount3D+1),jgT(mTerm3d),ggT(mTerm3d)
    if(node.gt.edgesNotTermCount3D) then
        i=node-edgesNotTermCount3D
        getValueOnEdge=0.d0
        do j=igT(i),igT(i+1)-1
            read(nFile,rec=jgT(j))v
            getValueOnEdge=getValueOnEdge+v*ggT(j)
        enddo
    else
        read(nFile,rec=node)getValueOnEdge
    endif
end
!*************************************************************************
subroutine basisFunctionOnVectPar(nEdge,x,y,z,hx,hy,hz,v,bfValue)
    dimension bfValue(3)
    real*8 bfValue,l0,l1,x,y,z,hx,hy,hz,v
    call basisFunctionOnElement(nEdge,x,y,z,bfValue)
    if(nEdge.ge.1.and.nEdge.le.4) then
        bfValue(1) = bfValue(1)*2.0/hx*v
    elseif(nEdge.ge.5.and.nEdge.le.8)then
        bfValue(2) = bfValue(2)*2.0/hy*v
    else
        bfValue(3) = bfValue(3)*2.0/hz*v
    endif
end
!*************************************************************************
real*8 function l0(x)
    real*8 x
    l0 = (1.d0 - x)/2.d0;
end
!*************************************************************************
real*8 function l1(x)
    real*8 x
    l1 = (x + 1.d0)/2.d0;
end
!*************************************************************************
subroutine basisFunctionOnElement(nEdge,x,y,z,bfValue)
    dimension bfValue(3)
    real*8 bfValue,l0,l1,x,y,z

    SELECT CASE (nEdge)
    CASE (1)
        bfValue(1) = l0(y)*l0(z)
        bfValue(2) = 0.d0
        bfValue(3) = 0.d0
    CASE (2)
        bfValue(1) = l1(y)*l0(z)
        bfValue(2) = 0.d0
        bfValue(3) = 0.d0
    CASE (3)
        bfValue(1) = l0(y)*l1(z)
        bfValue(2) = 0.d0
        bfValue(3) = 0.d0
    CASE (4)
        bfValue(1) = l1(y)*l1(z)
        bfValue(2) = 0.d0
        bfValue(3) = 0.d0
    CASE (5)
        bfValue(1) = 0.d0
        bfValue(2) = l0(z)*l0(x)
        bfValue(3) = 0.d0
    CASE (6)
        bfValue(1) = 0.d0
        bfValue(2) = l1(z)*l0(x)
        bfValue(3) = 0.d0
    CASE (7)
        bfValue(1) = 0.d0
        bfValue(2) = l0(z)*l1(x)
        bfValue(3) = 0.d0
    CASE (8)
        bfValue(1) = 0.d0
        bfValue(2) = l1(z)*l1(x)
        bfValue(3) = 0.d0
    CASE (9)
        bfValue(1) = 0.d0
        bfValue(2) = 0.d0
        bfValue(3) = l0(x)*l0(y)
    CASE (10)
        bfValue(1) = 0.d0
        bfValue(2) = 0.d0
        bfValue(3) = l1(x)*l0(y)
    CASE (11)
        bfValue(1) = 0.d0
        bfValue(2) = 0.d0
        bfValue(3) = l0(x)*l1(y)
    CASE (12)
        bfValue(1) = 0.d0
        bfValue(2) = 0.d0
        bfValue(3) = l1(x)*l1(y)
    END SELECT
end
!*************************************************************************
subroutine transformationCoordinats(x,y,z,hx,hy,hz,x1,x2,y1,y2,z1,z2,newX,newY,newZ)
    real*8 x,y,z,hx,hy,hz,x1,x2,y1,y2,z1,z2,newX,newY,newZ

    newX = (2.0*x - x1 - x2)/hx
    newY = (2.0*y - y1 - y2)/hy
    newZ = (2.0*z - z1 - z2)/hz
end
!*************************************************************************
!*************************************************************************!*************************************************************************